library("tidyverse")
library("here")
library("patchwork")

Prepare data

sensors <- read_csv(here::here("data/raw/sensors.csv")) %>% 
  filter(elevation == 1750)

d <- read_csv(here::here("data/raw/readings.csv")) %>% 
  filter(id_sensor %in% sensors$id_sensor)
## Warning: 794 parsing failures.
##  row            col expected actual                                                                                                  file
## 3418 positionsensor a double   NULL '/Users/ajpelu/Google Drive/MS/2021_MS_ECOSISTEMAS_MOUNTAIN/climate_qpyr_canar/data/raw/readings.csv'
## 3419 positionsensor a double   NULL '/Users/ajpelu/Google Drive/MS/2021_MS_ECOSISTEMAS_MOUNTAIN/climate_qpyr_canar/data/raw/readings.csv'
## 3468 positionsensor a double   NULL '/Users/ajpelu/Google Drive/MS/2021_MS_ECOSISTEMAS_MOUNTAIN/climate_qpyr_canar/data/raw/readings.csv'
## 3469 positionsensor a double   NULL '/Users/ajpelu/Google Drive/MS/2021_MS_ECOSISTEMAS_MOUNTAIN/climate_qpyr_canar/data/raw/readings.csv'
## 3470 positionsensor a double   NULL '/Users/ajpelu/Google Drive/MS/2021_MS_ECOSISTEMAS_MOUNTAIN/climate_qpyr_canar/data/raw/readings.csv'
## .... .............. ........ ...... .....................................................................................................
## See problems(...) for more details.
humedad <- d %>% 
  dplyr::filter(str_detect(tag, "Humedad")) %>% 
  dplyr::select(id_sensor, dat, time_dat, tag, positionsensor, tvariable) %>% 
  inner_join(sensors) %>% 
  dplyr::select(-tvariable, -elevation, -c(stationid:geo)) %>% 
  rename(pos = positionsensor) %>% as.data.frame()

Exploratory

h_dialy <- humedad %>% 
  group_by(date=lubridate::floor_date(time_dat, "day"), 
           id_sensor, tag, pos, location, replica) %>% 
  summarise(daily_avg = mean(dat)) %>% 
  filter(date > as.POSIXct("2017-12-31")) %>% 
  mutate(date = as.Date(date, format="%Y-%m-%d"))
           
colores <- c("CLARO" = "#CC7722", "ROBLEDAL" = "#0B6623")
ggplot(h_dialy, aes(x=date, y = daily_avg, colour = location)) + 
  geom_point(size=.5) + scale_color_manual(values=colores) + 
  facet_wrap(~tag) + theme_minimal()

h_dialy_avg <- h_dialy %>% 
  group_by(location, tag, date) %>% 
  summarise(daily_avg = mean(daily_avg))


ggplot(h_dialy_avg, aes(x=date, y = daily_avg, colour = location)) + 
  geom_point(size=.5) + scale_color_manual(values=colores) + 
  facet_wrap(~tag) + theme_minimal()

h_dialy_avg_spring <- h_dialy_avg %>% 
  mutate(m = lubridate::month(date)) %>% 
  filter(m %in% c(4:9))

ggplot(h_dialy_avg_spring, aes(x=date, y = daily_avg, colour = location)) + 
  geom_point(size=.5) + scale_color_manual(values=colores) + 
  facet_wrap(~tag) + theme_minimal() +
  scale_x_date(date_breaks = "3 months", date_labels = "%b", 
               sec.axis = dup_axis(labels = scales::date_format("%Y"), 
                                   breaks = scales::date_breaks("year")))

h_dialy_avg_spring %>% 
  filter(date > as.POSIXct("2018-12-31") &
           date < as.POSIXct("2019-12-31")) %>% 
  ggplot(aes(x=date, y = daily_avg, colour = location)) + 
  geom_point(size=.5) + scale_color_manual(values=colores) + 
  facet_wrap(~tag) + theme_minimal() +
  scale_x_date(date_breaks = "1 months", date_labels = "%b") +
  geom_line()

h_dialy_avg_spring %>% 
  filter(date > as.POSIXct("2017-12-31") &
           date < as.POSIXct("2018-12-31")) %>% 
  ggplot(aes(x=date, y = daily_avg, colour = location)) + 
  geom_point(size=.5) + scale_color_manual(values=colores) + 
  facet_wrap(~tag) + theme_minimal() +
  scale_x_date(date_breaks = "1 months", date_labels = "%b") +
  geom_line()

Datos Humedad suelo

  • Filtramos HumedadSuelo30profundidad
soil <- h_dialy %>% filter(tag == "HumedadSuelo30profundidad") 

ggplot(soil, aes(x=date, y = daily_avg, colour = location)) + 
  geom_point(size=.5) + scale_color_manual(values=colores) + theme_minimal()

  • Quitamos outliers (values > 50; values < 2)
soil <- h_dialy %>% filter(tag == "HumedadSuelo30profundidad") %>% 
  filter(daily_avg < 50 & daily_avg > 2) 

ggplot(soil, aes(x=date, y = daily_avg, colour = location)) + 
  geom_point(size=.5) + scale_color_manual(values=colores) + theme_minimal()

  • Valores agregados por replica
soil_avg <- soil %>% 
  group_by(date, location) %>% summarise(daily_avg = mean(daily_avg))
  • Valores de primavera y verano
soil_pv <- soil %>% 
  filter(lubridate::month(date) %in% c(3:9)) %>%
  mutate(jday = lubridate::yday(date)) %>% as.data.frame()

ggplot(soil_pv, aes(x=date, y = daily_avg, colour = location)) + 
  geom_point(size=.5) + scale_color_manual(values=colores) + theme_minimal()

  • Agrupamos por replica
# Ojo ver que pasa los días 170 y 164 en CLARO
# Parece que hay valores outliers en alguna replica 
soil_pv %>% filter(location == "CLARO") %>% 
  filter(jday %in% c(164, 170))
##          date id_sensor                       tag pos location replica
## 1  2018-06-13       117 HumedadSuelo30profundidad   1    CLARO      R1
## 2  2018-06-13       121 HumedadSuelo30profundidad   1    CLARO      R2
## 3  2018-06-13       123 HumedadSuelo30profundidad   1    CLARO      R3
## 4  2018-06-19       117 HumedadSuelo30profundidad   1    CLARO      R1
## 5  2018-06-19       121 HumedadSuelo30profundidad   1    CLARO      R2
## 6  2018-06-19       123 HumedadSuelo30profundidad   1    CLARO      R3
## 7  2019-06-13       121 HumedadSuelo30profundidad   1    CLARO      R2
## 8  2019-06-13       123 HumedadSuelo30profundidad   1    CLARO      R3
## 9  2019-06-19       117 HumedadSuelo30profundidad   1    CLARO      R1
## 10 2019-06-19       121 HumedadSuelo30profundidad   1    CLARO      R2
## 11 2019-06-19       123 HumedadSuelo30profundidad   1    CLARO      R3
## 12 2020-06-12       117 HumedadSuelo30profundidad   1    CLARO      R1
## 13 2020-06-18       117 HumedadSuelo30profundidad   1    CLARO      R1
##    daily_avg jday
## 1   23.00167  164
## 2   16.57000  164
## 3   13.22583  164
## 4   19.27583  170
## 5   12.92958  170
## 6   11.31125  170
## 7   12.30625  164
## 8    9.38600  164
## 9   13.58667  170
## 10  11.50375  170
## 11   7.34525  170
## 12  33.22000  164
## 13  45.87875  170
# Efectivamente valores anómalos en la R1 (en comparación con las otras réplicas). Los elimino 
soil_pv_avg <- soil_pv %>% 
  filter(!(location =="CLARO" & jday %in% c(164, 170) & replica == "R1")) %>% 
           group_by(date, jday, location) %>% summarise(daily_avg = mean(daily_avg)) 


ggplot(soil_pv_avg, aes(x=date, y = daily_avg, colour = location)) + 
  geom_point(size=.5) + scale_color_manual(values=colores) + theme_minimal()

  • Creamos perfil de esos meses
profile <- soil_pv_avg %>% 
  group_by(jday, location) %>% 
  summarise(mean = mean(daily_avg), 
            sd = sd(daily_avg), 
            se = sd/sqrt(length(daily_avg))
            ) %>% 
  filter(jday < 244) %>% 
  # ojo esta fecha es simplemente para visualizacion
  mutate(date = as.Date("2018-01-01") + jday -1)

p_profile <- ggplot(profile, aes(x=date, y = mean, colour = location, fill=location)) + 
  # geom_point(size=.5) + 
  geom_line() + 
  geom_ribbon(aes(ymin=mean-se, ymax=mean+se), colour=NA, alpha =.4) + 
   # geom_errorbar(aes(ymin=mean-se, ymax=mean+se))
  scale_color_manual(values=colores, labels = c("Open","Forest")) +
  scale_fill_manual(values=colores, labels = c("Open","Forest")) +  
  theme_bw() + 
  theme(legend.position = "bottom",
        legend.title = element_blank(), 
        panel.grid = element_blank()) +
  scale_x_date(date_breaks = "1 month", date_labels = "%b") + 
  xlab("") + 
  ylab("Contenido Volumétrico Agua (-30 cm) (%)")

p_profile

ggsave(plot = p_profile,
       filename = here::here("figs/profile_humedad_verano.jpg"), 
       height = 5, width = 5)
profile1819 <- soil_pv_avg %>% 
  filter(lubridate::year(date) %in% c(2018:2019)) %>% 
  group_by(jday, location) %>% 
  summarise(mean = mean(daily_avg), 
            sd = sd(daily_avg), 
            se = sd/sqrt(length(daily_avg))
            ) %>% 
  filter(jday < 244) %>% 
  # ojo esta fecha es simplemente para visualizacion
  mutate(date = as.Date("2018-01-01") + jday -1)


ggplot(profile1819, aes(x=date, y = mean, colour = location, fill=location)) + 
  # geom_point(size=.5) + 
  geom_line() + 
  geom_ribbon(aes(ymin=mean-se, ymax=mean+se), colour=NA, alpha =.4) + 
   # geom_errorbar(aes(ymin=mean-se, ymax=mean+se))
  scale_color_manual(values=colores) +
  scale_fill_manual(values=colores) + 
  theme_minimal() +
  scale_x_date(date_breaks = "1 month", date_labels = "%b") + xlab("") + ylab("Humedad Suelo (-30 cm) (%)")

Datos de Precipitación para esos tres años (2018 a 2020)

  • Datos de prec solo hasta enero de 2020
dfprec <- read_csv(here::here("data/raw/psn07_prec.csv")) %>% 
  mutate(date = as.Date(date, format = "%Y-%m-%d"))

prec <- dfprec %>% filter(lubridate::year(date) %in% c(2018:2019)) 

p <- ggplot(prec, aes(x=date, y = value)) + 
  # geom_density(stat = "identity", fill="blue", colour = "blue", 
  #              outline.type = "upper", 
  #              size = .8) + 
  geom_bar(stat = "identity", colour=NA, fill= "blue", width = 3) + 
  theme_bw() + 
  theme(
    panel.grid.minor =element_blank(),
    axis.title.y =  element_text(size = 8)
  ) + 
  ylab("Precipitación (mm)") + xlab("") +
  scale_x_date(position = "top")
  • combinar
s <- soil_avg %>% 
  filter(lubridate::year(date) %in% c(2018:2019)) %>% 
  ggplot(aes(x=date, y = daily_avg, colour = location)) + 
  # geom_point(size=.5) + 
  scale_color_manual(values=colores,
                     labels = c("Open","Forest")) + 
  geom_line() + 
  theme_bw() +
  theme(
    panel.grid.minor =element_blank(), 
    legend.position = "bottom", 
    legend.title = element_blank(),
    axis.title.y =  element_text(size = 8)
  ) + 
  ylab("Contenido Volumétrico Agua (%)") + xlab("")


plot_conjunto <- p/s 
plot_conjunto
## Warning: position_stack requires non-overlapping x intervals

ggsave(plot = plot_conjunto,
       filename = here::here("figs/profile_humedad.jpg"),
 height = 5, width = 6)
## Warning: position_stack requires non-overlapping x intervals